######################################################################
#########            FEMALE EMPOWERMENT INDEX OVER TIME   ############
######################################################################
load("Figuredata.Rdata")

#1900-2015
mydata2 <- subset(mydata, year>1899 & mydata$year<2016)
avg <- with(mydata2, aggregate(v2x_gender, by=list(mydata2$year), FUN=mean, na.rm=T))
avg2 <- with(mydata2, aggregate(v2x_gender, by=list(mydata2$year), FUN=sd, na.rm=T))

pdf("fig1.pdf", width=15, height=10)
plot(mydata2$year, mydata2$v2x_gender, type="n",
     xlab="Year", ylab="Avg. global female empowerment index (V-DEM)",
     col="grey")
lengths <- seq(from=1900, to=2015, by=1)

avg2$sd1 <- avg$x+avg2$x
avg2$sd2 <- avg$x-avg2$x

polygon(c(lengths
          , rev(lengths))
        , c(avg2$sd1
            , rev(avg2$sd2)),
        col="lightblue", density=30,border=NA)
lines(avg$Group.1, avg$x, col="black", lwd=3)
dev.off()


#1817-2015
mydata <- subset(mydata, year<1900)
avg <- with(mydata, aggregate(hist_v2x_gender, by=list(mydata$year), FUN=mean, na.rm=T))
avg2 <- with(mydata, aggregate(hist_v2x_gender, by=list(mydata$year), FUN=sd, na.rm=T))

pdf("fig2.pdf", width=15, height=10)
plot(mydata$year, mydata$hist_v2x_gender, type="n",
     xlab="Year", ylab="Avg. global female empowerment index",
     col="grey")
lengths <- seq(from=1790, to=1899, by=1)

avg2$sd1 <- avg$x+avg2$x
avg2$sd2 <- avg$x-avg2$x

polygon(c(lengths
          , rev(lengths))
        , c(avg2$sd1
            , rev(avg2$sd2)),
        col="lightblue", density=30,border=NA)
lines(avg$Group.1, avg$x, col="black", lwd=3)
dev.off()


